Modeling and probabilistic analysis of civil aircraft operational risk for suborbital disintegration accidents

To reduce the collision risk to civil airliners caused by suborbital vehicle disintegration events, this paper uses a covariance propagation algorithm to model the debris landing point of suborbital disintegration accidents and gives a collision probability analysis method for civil airliners encountering debris during the cruise. Collision warning is performed for airborne risk targets to improve the emergency response capability of the ATC surveillance system to hazardous situations. The algorithm models the three-dimensional spatial motion target localization problem as a Gauss-Markov process, quantifying the location of debris landing points in the vicinity of nominal trajectories. By predicting the aircraft trajectory, the calculation of the inter-target collision probability is converted into an integration problem of a two-dimensional normally distributed probability density function in a circular domain. Compared with the traditional Monte Carlo method, the calculation speed of debris drop points is improved, which can meet the requirements of civil aviation for real-time response to unexpected situations.


Introduction
In recent years, along with the shift from government-driven to market-driven drivers of space activities, commercial spaceflight has become an important driver of global space activities, and the commercialization of suborbital flight has developed rapidly [1]. Unlike civil aviation flight methods, the suborbital flight is a flight between both gravity and aerodynamics, as the speed of flight around the earth is not reached, the vehicle is equivalent to making a parabola back to earth in space [2]. Compared to a normal flight, the suborbital flight time will be shortened several times; compared to a spacecraft, suborbital vehicles are cheap and reusable, thus presenting a huge economic value [3]. It is predicted that suborbital transportation is likely to become one of the main modes of space travel in the next generation [4].
Suborbital vehicles travel at hypersonic speeds, with orbits between the maximum altitude of existing aircraft and the minimum orbital altitude of satellites, and are prone to disintegration and large debris generation under strong aerodynamic loads. The Inter-Agency Space Debris Coordination Committee (IADC) Space Debris Mitigation Guidelines [5]  spacecraft typically disintegrate at an altitude of 75 km, and that a typical spacecraft disintegration altitude of (80±10) km has been observed in practice, with object re-entry survival The aerodynamic forces on these fragments in the atmosphere are somewhat random and their distribution and landing points are difficult to predict, with serious consequences in the event of a collision with a civil aircraft [6]. During space launches, traffic control authorities completely close off large, infinitely high areas of "hazardous airspace [7]" in advance of the launch plan. This measure avoids aircraft encounters with risky objects and maximizes the safety of civil aviation operations, but also causes extensive flight delays and economic losses. As the number of commercial suborbital space activities increases in the future, this static, excessive, and prolonged approach to traditional airspace management will be costly to airlines [8]. Suborbital debris distribution modeling is a key technique to form the danger zone prediction, and effective danger zone prediction is a prerequisite for implementing safe air traffic control. In terms of debris distribution modeling, some scholars have approximated the suborbital debris prediction hazard zone as a rectangle based on the debris statistics of the Columbia disintegration event [9][10][11], with the length of the rectangle being the disintegration height divided by 1000 and the width of the rectangle being 1/8 of the length; this approximation is relatively crude, lacks rigorous scientific justification, and does not have universal significance. Some scholars have used the TAP (Trajectory Analysis Program) model to simulate the debris distribution in a visual simulation system for suborbital accident debris risk assessment [12,13], which was first designed to simulate the debris distribution generated by the accidental disintegration of an aircraft in an aerobatic stunt show and is not essentially applicable to the trajectory prediction of suborbital disintegration accident debris. In addition, in risky target collision warnings, some scholars use the orbital state and error covariance information obtained from the forecast to calculate the collision probability between targets, but this algorithm is only applicable to outer space targets.
To improve the real-time response capability of ATMS to suborbital disintegration accidents and minimize the waste of airspace resources, this paper uses a covariance propagation algorithm to model the suborbital disintegration accident debris position distribution as a Gauss-Markov process, combined with the classical aircraft track prediction model, fully considering the safety requirements in the civil aviation field and the response speed requirements of ATC system, the collision between aircraft and suborbital debris The probability is calculated, and the calculation results can provide real-time reference information for air traffic controllers to issue redirect instructions to pilots for air risk target collision warning.

Covariance propagation method for debris model
Reyhanoglu et al. [14] proposed a covariance propagation algorithm for determining the fourdimensional footprint of debris based on the three-degree-of-freedom (3-DOF) model of falling objects on a rotating planet. The algorithm is based on aerodynamics and Newton's laws of motion to establish a mathematical model of the diffusive motion of individual debris in a complex atmospheric environment, and further analysis is performed. Let x 1 x 2 x 3 denote the position of the debris in the station center coordinate system at the moment of aircraft disintegration. The origin of the coordinate system is at (θ 0 , ϕ 0 ) the surface of the earth, where (θ 0 , ϕ 0 ) is the initial latitude and longitude at the moment of aircraft disintegration. The x 1 axis points to the east and the x 2 axis points to the north, the x 1 x 2 plane represents the plane tangent to the earth at the origin, and the x 3 axis is perpendicular to the plane pointing to the zenith direction. The coordinate system can also be called the ENU (East-North-Up) coordinate system. Then the equation of motion describing the propagation trajectory of the debris can be given by the following equation: Where: 3 ] T represent the position vector and velocity vector in the ENU coordinate system, respectively; ω = [0,ω e cos ϕ 0 , ω e sin ϕ 0 ] T is the rotational angular velocity vector in the ENU coordinate system, where ω e = 7.2921 × 10 −5 rad/s represents the mean angular velocity of the Earth's rotation; e 3 = [0, 0, 1] T represents the third standard unit vector; ξ represents the random acceleration vector caused by modeling uncertainty and disturbances; R e = 6378 × 10 3 m represents the radius of the Earth, and according to the inverse square gravity model, the gravitational acceleration can be represented by g = g e (R e /(R e + x 3 )) 2 , where g e = 9.81/s 2 ; a D represents the instantaneous acceleration associated with the atmospheric density, which can be expressed as a D t , where β represents the ballistic coefficient, and v a is the difference between the velocity vector v of the debris and the wind vector w(x,t), v a = vw(x, t).a G represents the acceleration of gravity, a G ¼ gê 3 .The atmospheric index model can be expressed as r x 3 Þ ¼r e exp À x 3 H À � À . For debris generated from suborbital disintegration accidents, the Monte Carlo method is usually used to simulate the Spatio-temporal evolution process of their state vectors. However, in practical applications, the simulation process of a large amount of debris is time-consuming and cannot meet the high requirements of real-time response for air traffic control work. Therefore, the estimation of debris drop-off points using the Covariance Analysis Description Equation Technique (CADET) can significantly reduce the computational effort. The algorithm describes the stochastic process of predicting debris drop points around a linearized nominal trajectory as a Gauss-Markov process and then constructs a probability ellipsoid of the Gauss-Markov process at a certain confidence level through a probability density function to quantify the location of debris drop points near the nominal trajectory.
First, define the extended state vector s = [x T , v T ] T , then the fragment motion Eq (1) can be reformulated as: Where: η = [0 T , ξ T ] T . Define _ s � as the nominal state vector, i.e., _ s � ¼ f ðs � Þ, and the perturbation vector z = s-s�, then Eq (2) is linearized as: Where Assuming that the initial state obeys the Gaussian distribution and that the noise statistics are characterized by

Þ, then Eq (3) can be regarded as a continuous
Gauss-Markov process. The fragment position vector at a time t can be expressed as a Gaussian random variable with mean x t In the above model, the main uncertainty influences include ballistic coefficients, atmospheric density, and wind vectors. The widely used empirical density model MSISE-00 (Mass Spectrometer Incoherent Scatter Radar Extended 2000 [15]) and the wind model HWM-07 (Horizontal Wind Model 2007 [16]) were chosen for modeling the disintegration accident. The debris ballistic coefficients can be characterized by probabilistic statistics. The motion state (longitude, latitude, altitude, velocity vector, heading angle, etc.) of the suborbital vehicle before disintegration can usually be derived from radar observations or flight trajectory equations. The description of the debris dispersion characteristics is extended to four dimensions by introducing an ellipsoidal set constructed by parameterizing the center and shape matrices concerning time. The resulting 4D probability constrained problem is much more complex than the 3D problem due to the infinite nature of the time variables. Therefore, solving the problem requires discretizing the time and associating each sample time with an ellipsoid. In this case, the 4D footprint of a single momentary fragment is part of the probability ellipsoid, while a series of instantaneous probability ellipsoids form a "fragment channel", as in Fig 1.

Aircraft track prediction model
The coordinate system describing the position of the moving target mainly includes the carrier coordinate system and the navigation coordinate system, and the target attitude is described by the pitch angle θ, roll angle ϕ, and yaw angle ψ. The carrier coordinate system is the Front-Left-Up (FLU) reference system, and the navigation coordinate system is the ENU reference system, as shown in Fig 2, and the angular conversion relationship between the carrier coordinate system and the navigation coordinate system is represented by the Earth Centered Earth Fixed (ECEF) reference system.
In this paper, we only discuss the case that the aircraft encounters suborbital debris during the cruise (the cruise altitude of civil aircraft is generally around 10km), i.e., the pitch angle is considered to be 0 and the flight attitude is neglected. As shown in Fig 3, let be the angle between the aircraft heading and due east direction, and let R n b denote the coordinate conversion matrix from the carrier coordinate system to the navigation coordinate system. Due to the presence of uncertainties such as wind speed, navigation, and aircraft positioning errors, after studying and analyzing a large amount of track data, it is generally assumed that the track prediction error obeys normal distribution. Paielli et al. [17] proposed that the position error between the aircraft's actual position within the horizontal altitude layer and the planned track within the next 30 minutes follows a zero-mean normal distribution and that the position prediction errors are independent of each other along the track direction and in the vertical track direction. In an aircraft equipped with Flight Management System (FMS), the aircraft adjusts the navigation position through a feedback mechanism. In the vertical track direction, the position closed-loop feedback can stabilize the Root Mean Square (RMS) error at a constant level, while in the along-track direction, due to the unpredictable wind direction, the FMS can only achieve feedback through airspeed, resulting in a classical linear growth of the along-track RMS error [18], and the along-track error and the vertical track error are independent of each other of the aircraft, the aircraft track error can be obtained as an ellipse with its long axis along the flight direction and its short axis perpendicular to the flight direction, as shown in Fig 3. The position prediction error of the aircraft in the carrier coordinate system is a diagonal matrix. Suppose the aircraft position in the carrier coordinate system is q, and the corresponding predicted position is � q, and the trajectory prediction errorq ¼ q À � q. The track prediction error approximately conforms to the normal distribution, and its covariance matrix is the diagonal matrix S.
Where: σ x is the error along the heading direction, when the aircraft travels at a uniform speed in a straight line, σ x increases linearly with time, and satisfies: Where: r a = 0.25n.mi/min is the growth rate of error along with the heading.σ y is the error along the vertical heading direction and satisfies: Where: r c = 1/57 represents the error growth factor and s(t) represents the distance flown by aircraft in the predicted time. σ z is the error along the vertical heading direction and is usually assumed to be σ z = 30 m.

Definition of the encounter coordinate system
To avoid the impact of vehicle debris in suborbital disintegration accidents on the safe operation of civil airliners, the collision probability of debris and aircraft in the same space needs to be considered comprehensively for airborne risk target collision warning. To simplify the calculation process of collision probability, the calculation is generally based on the following assumptions [19]: the position vector and velocity vector of two targets in the navigation coordinate system during the encounter are known and can be equivalently transformed into a sphere of known radius; the motion process of two targets during the encounter has the characteristics of uniform linearity and no uncertainty in velocity, i.e., the position error ellipsoid is guaranteed to remain constant during the encounter; the position error of two targets obeys the normal distribution in three-dimensional space, which can be described by the predicted position and position error covariance matrix. The collision condition is that the distance between the two targets is less than the sum of their equivalent radii, and thus the collision probability is defined as the probability that the mode of the relative position vectors is less than the sum of their equivalent radii. The position-velocity geometric relationship between the aircraft and the individual debris after a suborbital disintegration accident is shown in Fig 4. Where v 1 , v 2 denotes the velocity vectors of the two targets, respectively, r 1o , r 2o is the predicted position vector, and r 1 , r 2 is the actual position vector. The position errors of the two targets are expressed in the carrier coordinate system as e 1 , e 2 , and the equivalent radii are R 1 , R 2 , respectively. The actual position vector of the target can be expressed as the predicted position vector plus the position error vector, i.e., r 1 = r 1o + e 1 , r 2 = r 2o + e 2 .
Let P c denote the collision probability and the coordinate transformation matrix of the navigation and carrier coordinate systems is known. The distance between two targets is the relative position vector modulo ρ = |ρ| = |r 1 -r 2 |. The collision probability is expressed as the probability that the distance between two targets is less than the sum of equivalent radii, so the collision probability P c can be expressed as P c = P {ρ � R = R 1 + R 2 }. Let the starting moment be t 0 = 0, then the relative position vector between the two targets at the moment t is: Where: v r is the relative velocity vector. The square of the distance between the two targets is ρ 2 (t) = ρ(t) � ρ(t), find the derivative concerning time and let the derivative be zero as follows: From Eq (9), the minimum distance between the two targets, i.e., the time to reach the Closest Point of Approach (CPA), is: Then, the relative position vector of the two objects is: Further: That is, when the distance between two targets is smallest, the dot product of the relative position vector ρ(t cpa ) and the relative velocity vector v r is 0, which means that the two vectors are perpendicular to each other. In other words, when two objects are closest, they are in the plane perpendicular to the relative velocity vector, which is defined as the Encounter Plane (EP). The encounter coordinate system o-x e y e z e is defined as shown in Fig 5, where the origin o 2 of the coordinate system is located at the distribution center of the aircraft (target 2), the x e axis and y e axis are within the EP and are perpendicular to each other, the x e axis points to the projection point of the distribution center of the debris (target 1) on the EP, and the z e axis point to the direction of the relative velocity vector. According to the definition of the encounter coordinate system, the unit vectors of the three axes of the encounter coordinate system are: Thus, the coordinate transformation matrix M e for the navigation coordinate system to the encounter coordinate system:

Position error projection
Suppose a vector is represented in the carrier coordinate system, the navigation coordinate system, and the EP coordinate system as r b , r n , r e , respectively, then the transformation relation is: The resulting coordinate transformation matrix from the carrier coordinate system to the encounter coordinate system is M e1 ¼ M e R n b . Let X r denote the random vector in the carrier coordinate system, which is represented in the EP coordinate system X re , and the transformation relationship between them can be expressed as X re = M e1 X r . Then the variance matrix X re in the encounter coordinate system is: This converts the position error covariance matrix represented in the carrier coordinate system to the position error covariance matrix represented in the encounter coordinate system, and by eliminating the terms associated with z e in this covariance array, the position error is projected onto the encounter plane, thus reducing the three-dimensional complex problem to a two-dimensional one, as shown in Fig 6: The position vectors X 1 and X 2 of two targets in the encounter plane obey a two-dimensional normal distribution and are independent of each other, then the relative position vector X = X 1 -X 2 is also a normal random vector. The collision probability of two targets is the probability that the relative position vectors fall into the circle with R = R 1 + R 2 as the radius. Therefore, the joint circular domain is constructed with the origin o 2 as the center and R as the radius, and the position error ellipses of the two targets are formed into a joint error ellipse at target 1, as shown in Fig 7: In general, the variance matrix of the relative position vector X is not diagonal, i.e., the EP coordinate axes are not parallel to the direction of the principal axes of the joint error ellipse, and to simplify the calculation, the coordinate transformation is needed again. To make the principal axes of the joint error ellipse parallel to the coordinate axes, define a computational coordinate system o-x c y c , with the origin coinciding with the encounter coordinate system, the x c axis pointing in the direction of the short axis of the joint error ellipse in the encounter plane, and y c pointing in the direction of the long axis. Assuming that an angle of rotation θ is required, in the encounter co-ordinate system, let the variance array X be where: κ denotes the correlation coefficient. In the computational coordinate system o-x c y c , the variance matrix of the random vector X is Based on the relationship between the variance array under coordinate transformation, it follows that Where: M(θ) is the coordinate transfer matrix from the encounter coordinate system to the calculated coordinate system. Since Var (X) c is a diagonal array, it follows that For the axis of rotation x c to point in the direction of the short semi-axis of the error ellipse, there should be s The angle of rotation θ should therefore satisfy the condition that The distribution center should rotate in the computational coordinate system by an amplitude of θ' = -θ, as shown in Fig 8, at which point the parameters of the error distribution in the computational coordinate system are The two-dimensional normal distribution probability density function (PDF) is: The collision probability P c is the integral of the PDF in the joint circular domain: In this way, the problem of calculating probability is transformed into the problem of integrating the probability density function in the circle domain.

Example analysis
Statistical data show that the aircraft in cruise is in a relatively stable state with typical motion characteristics, and its motion is considered to satisfy a typical linear growth pattern without additional assumptions, where the typical airspeed of a commercial transport aircraft is 0.25 km/s. The motion of the vehicle before the disintegration accident is usually known from radar observations and flight trajectory equations, and the international common spacecraft The disintegration altitude is 78 km, and the yaw angle ψ 0 = 0, pitch angle θ 0 = -1, initial velocity v0 = 7.1km/s, and initial longitude and latitude of the suborbital vehicle before disintegration are assumed to be γ 0 = 157˚W and φ 0 = 20˚N, respectively. The aircraft model takes the Boeing 737 as an example, and the aircraft's equivalent radius can be set to 30 m. In this modeling, only the individual debris form is considered, and the suborbital debris equivalent radius is set to 0.2 m based on the debris statistics of the Columbia disintegration accident, and the joint equivalent radius is R = 30.2 m. Through the above calculation, when the debris of the suborbital disintegration accident falls to the civil aviation altitude, the state information of the aircraft and some debris at the reference moment is given in combination with the aircraft trajectory prediction model, and the collision probability, the closest moment and the closest distance between the aircraft and each debris are calculated respectively, as shown in Table 1.

Debris position error distribution
It is assumed that the error of each debris position is spherically distributed, that is: Take the first set of simulation data and calculate the collision probability between the aircraft and debris 01, as shown in Fig 9. The collision probability P c increases first with the increase of position error and starts to decrease after σ x = 103 m reaching a great value P cmax . When the position error is small, as long as the aircraft and the debris are not in the intersecting path, the probability of collision is small. As the position error increases, the probability of collision increases, and after reaching the maximum value, as the error range continues to increase, the probability of collision is small even if the distance between the aircraft and the debris is small. In addition, in the actual situation the debris position error is larger in the direction of the z axis, therefore, considering the position error distribution σ x : σ y : σ z increases from 1:1:1 to 1:1:5, as seen in Figs 10 and 11, the collision probability extreme value point keeps shrinking.

Airborne risk target collision warning
The international general spacecraft maneuver avoidance probability yellow threshold value is p c = 10 −5 , the red threshold is p c = 10 −4 , according to the calculation results in Table 1, it is known that the collision probability between the aircraft and debris 01 and 02 are greater than and equal to the red threshold value, respectively, and the collision probability with debris 03 is in the yellow warning level, there is a possibility of collision and need to take avoidance maneuvers in advance. In addition, the collision probability between the aircraft and debris 04 is much smaller than the yellow warning value and will be a safe rendezvous event. In summary, to avoid the safety impact of suborbital disintegration on the civil airliner, the collision probability value between the aircraft and any debris in the debris group should be less than the yellow warning value, and the collision probability between the aircraft and the debris should be reduced to below 10 −7 after the avoidance measures are taken.

Fragmentation prediction calculation efficiency
The covariance propagation method used in this paper has a more significant advantage in running time compared to the traditional Monte Carlo algorithm while increasing the sampling points, as shown in Fig 12, which will facilitate the real-time response work of ATMs.

Conclusion
In this paper, the covariance propagation algorithm is used to model and analyze the debris drop points generated by suborbital vehicle disintegration accidents. Compared with the traditional Monte Carlo algorithm, the computational efficiency of this algorithm is significantly improved. Combined with the aircraft trajectory prediction model, the position and velocity information of the aircraft and each suborbital debris at the closest moment in the same altitude level are predicted in advance, and the position error covariance matrix of the aircraft and suborbital debris are combined by the above method. After linear transformation and coordinate conversion, the problem of collision probability calculation is converted into the problem of integrating the probability density function of two-dimensional normal distribution in a circular domain. Concerning the international common standard for spacecraft avoidance maneuvers, the collision warning of airborne risk targets is carried out, and an important foundation for the construction of a safe and efficient decision support system for flight avoidance is laid.